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Abstract 

I consider two coupled oscillators, each modulating 
the other’s frequency. This system is governed by 
four parameters: the base frequency and modulation 
index for each oscillator. For some parameter values 
the system becomes unstable. I use the Lyapunov 
exponent to measure the instability. I generate im¬ 
ages of the parameter space, implementing the num¬ 
ber crunching on graphics hardware using OpenGL. 
I link the mouse position over the displayed image 
to realtime audio output, creating an audio-visual 
browser for the 4D parameter space. 
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1 Introduction 

In Soft Rock EP [ClaudiusMaximus, 2005] and 
Soft Rock DVD [ClaudiusMaximus, 2006] I ex¬ 
plored the transitions between order and chaos 
in coupled FM oscillators, implemented in Pure- 
data using GridFlow for visualization of the out¬ 
put waveforms over time. More recently I devel¬ 
oped the idea to map the parameter space on a 
perceptually relevant level and in performance 
choose parameters on the basis of desired sound 
character. 

Seeing a bifurcation diagram produced by an 
analogue Moog synthesizer [Slater, 1998] and 
images of Lyapunov fractals [Dewdney, 1991], 
I decided to apply the latter technique to cou¬ 
pled FM oscillators in the digital realm. 



Figure 1: Example output. 



I had originally experimented with one Pure- 
data batch mode instance per CPU core each 
sending analysis data to a realtime Pure-data 
instance. The analysis used various methods 
(including FFT for spectral statistics and the 
sigmund external for pitch tracking) to clas¬ 
sify points into pitched (ordered, stable) or un¬ 
pitched (chaotic, unstable) with measures of 
distortion or noisiness. Sadly this approach 
proved impractical as it achieved only tens 
of pixels per second, even with a fast multi¬ 
core CPU, and porting these signal analysis 
algorithms to massively-parallel programmable 
graphics hardware seemed to be too difficult. 

2 Formulation 

2.1 Coupled FM Oscillators 
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Figure 2: Coupled FM oscillators in Pure-data. 


I’ll write the four-dimensional parameter 
space vector 
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and the (2d+2)-dimensional phase space vector 
Z = (%ni Uni %n— lj Un— 1 1 • • • > 3-n—di Vn—d) 

with sample rate SR = 48000. I’ll fix d = 1 for 
reasons explained in Section 5.1. 

2.2 Lyapunov Exponents 

A good introduction is found in Chapter 4.3 
Lyapunov Exponent [Elert, 2007]. The defini¬ 
tion is covered in Chapter 13.7 Liapounov expo¬ 
nents and entropies [Falconer, 2003] which also 
relates it to measures of fractal dimension. 

The Lyapunov exponent A measures diver¬ 
gence in phase space: 
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An attracting orbit has A < 0 and a divergent 
(chaotic) orbit has A > 0. 

The wrapping of phase into [0,1) requires a 
modified norm: 
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Consider the two coupled oscillators in Figure 2. 
I formulate this into a mutual recurrence rela¬ 
tion: 

x n +i = %(x n + I(/ x + m x cos(27r y n - d ))) 

Un+l = %{y n + I (fy + m y COs(27 TX n - d ))) 

where 
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Here x n , y n is the phase of each oscillator at 
time step n, d is a delay measured in samples, 
f x , f y is the base frequency of each oscillator as 
a MIDI note number, and m x , m y is the mod¬ 
ulation index of each oscillator as a MIDI note 
number. 


2.3 Viewing Planes 

An image is 2D, which requires choosing a sub¬ 
set of the 4D parameter space to visualize. I 
chose two particular planes: 


A + (a 0 ,r 0 ) 


A-(a 0 ,r 0 ) 


( 


do + A) 


V 

( 


ao + A) 


V 


1 

1 

0 

0 

1 

-1 

0 

0 


0 

0 

1 

1 

0 

0 

1 

-1 


u 

v 


u 

v 


( 3 ) 

where (it, v) is the coordinates of the pixel, ao 
is the centre of the view, and ro is the radius of 
the view. 





3 Results 

3.1 Examples 

Figure 3(a) shows the initial view on starting 
the interactive browser. Low frequencies to the 
left are stable even at high modulation index 
away from the central axis. High frequencies to 
the right become chaotic at progressively lower 
modulation index, (b) shows the plane at 
the same location, (c) shows bands alternat¬ 
ing between stability and chaos. The bands be¬ 
come distorted and collapse as the modulation 
index and frequency increase, (d) shows its 


plane, bands become rings. When the frequency 
is greatly increased, the shapes become more 
intricate, (e) exhibits spirals of stability, with 
similar spirals in the plane in (f). 

When f x = f y and m x = m y the A + plane 
has mirror symmetry about its horizontal axis, 
and the plane has two-fold rotational sym¬ 
metry about its centre. Breaking the symmetry 
and setting f x ^ f y or m x / m y leads to di¬ 
verse forms, shown in Figure 4. In particular 
Figure 4(c) has shapes that resemble those of 
Lyapunov space images of the logistic map. 






(a) A_((117.0,148.4, 20.4, 2.7), 1.8) 



(b) A_ ((141.46,146.22, 22.76, 0.27), 0.14) 



(c) A+((103.65,108.41,33.42,10.93),0.14) 



(d) A- ((89.8,137.5, -17.5, -7.1), 3.7) 
Figure 4: More examples. 


3.2 Interactive Explorer 

I implemented an interactive audio-visual ex¬ 
plorer for the parameter space of coupled FM 
oscillators. Clicking with the mouse zooms the 
view about the clicked point. The left button 
(or scroll up) zooms in, the right button (or 
scroll down) zooms out, the middle button cen¬ 
ters the view on the target point. Pressing the 
TAB key toggles between the A+ and A- planes 
in Equation 3, and Fll toggles full screen oper¬ 
ation. 

While the GPU simulates and analyses one 
oscillator pair per pixel, the CPU simulates one 
oscillator pair with a determined from the pixel 
under the mouse pointer. The image acts as a 
map, a reference frame for chosing parameters 
to audition by moving the mouse. 

4 Implementation 

I used OpenGL [Segal, 2013] and OpenGL 
Shading Language [Kessenich, 2013], with 
GLUT [Kilgard, 1996] for windowing and input 
event handling, and JACK [Davis, 2013] for au¬ 
dio output. 

4.1 Introduction to Modern OpenGL 

Modern OpenGL has a programmable shader 
pipeline. Vertex attributes are read from vertex 
buffers and processed by vertex shaders. The 
outputs of the vertex shader (called varyings) 
are further manipulated by an optional geome¬ 
try shader stage. Geometry shaders can output 
a different vertex count to their input count, 
whereas vertex shaders are one-in one-out. The 
result of the geometry shader can be captured 
into another vertex buffer using transform feed¬ 
back. Following the geometry shader the prim¬ 
itives (points or triangles) are rasterized, and 
varyings interpolated across each primitive. Fi¬ 
nally a fragment shader takes these values and 
computes the colour at that pixel. The output 
of a fragment shader can be captured by attach¬ 
ing a texture to a framebuffer. 

4.2 Computation Overview 

To render an image I first fill a texture with 
(u, v) coordinates using a framebuffer object 
and a fragment shader. I copy this texture 
to a vertex buffer object, interleaving an initial 
phase space vector z = (0,0,0,0) and Lyapunov 
exponent statistics vector l = (0,0, 0,0). 

Using a vertex shader, I calculate a from 
(u, v) using Equation 3, and then compute a 
rough estimate of the Lyapunov exponent us¬ 
ing Equation 2 by perturbing zi(0) = zq(0) + 5 


with 5 small and performing t = 256 iterations 
of Equation 1. I discard the first few repetitions 
and those resulting in — oo, and accumulate the 
rough A estimates in l. 

Between each repetition I compact the work¬ 
ing set using a geometry shader, plotting 
points whose mean Lyapunov exponent esti¬ 
mate changed very little during the previous 
step. I keep the others to refine further, direct¬ 
ing the computational effort on the points that 
need it most: those slow to converge. 

To ensure user interface responsiveness, the 
computation is amortized over several frames. 
I divide the target frame period by the mea¬ 
sured time for one repetition to compute how 
many repetitions to perform that frame. The 
repetitions-per-frame increases as the working 
set becomes smaller. 

4.3 Noise Increases Stability 

At the end of each repetition I keep z\ instead 
of zq. This effectively adds a small amount 
of noise, counter-intuitively increasing stability. 
Noise allows more of the phase space to be ex¬ 
plored, and makes it more likely for the per¬ 
turbed orbit to reach an attracting part of the 
phase space. 

4.4 Dither Increases Quality 

To reduce grid sampling artifacts, I perturb 
(u, v) within the bounds of its corresponding 
pixel before calculating the a parameter vector 
for each repetition. 

5 Conclusions 

5.1 OpenGL Issues 

The current implementation is hardcoded with 
delay d = 1 and would be very awkward to 
generalize. OpenGL architecture limits each 
vertex attribute to four components with the 
maximum number of attributes typically lim¬ 
ited to sixteen. This totals 64 floats per ver¬ 
tex, 6 of which are needed for the pixel coor¬ 
dinates and Lyapunov exponent statistics accu¬ 
mulation. Therefore using OpenGL imposes a 
limit d < 28. For comparison my original ex¬ 
periments in Soft Rock EP used Pure-data’s de¬ 
fault block size of 64, with d = 32. Moreover, 
increasing d increases video memory consump¬ 
tion. With the maximum d = 27, browsing at 
1920 x 1080 resolution would require over 1GB. 

My future work on this project will look into 
using OpenCL, which provides a heterogenous 
CPU and GPU computation framework. I hope 


it will avoid the inherent awkwardness of abus¬ 
ing OpenGL shaders to perform calculations. 

5.2 Audio Issues 

While the implementation works as intended, 
with d = 1 the sound is nowhere near as rich 
and varied as with d = 32. With small d there 
is much more very high frequency content in 
interesting-looking regions. I’ve not found any 
regions of the parameter space with both in¬ 
teresting appearance and palatable audio fre¬ 
quencies at d = 1, while with high d there are 
parameters that generate sounds that fluctuate 
intermittently between smooth tones and noise. 

I haven’t been able to visualize with high d to 
know whether their neighbourhoods look as in¬ 
teresting as they sound. 

Unfortunately, heavy use of the GPU in the 
interactive browser can block the operating sys¬ 
tem for too long and cause audible glitches 
(JACK xruns). This situation may change as 
free drivers continue to improve, allowing use of 
the browser in a live situation. 

5.3 Pretty Pictures 

Despite these shortcomings, I think the images 
look good. I plan to render a selection at high 
resolution and print postcards and posters. For 
huge images I divide the image plane into tiles 
and compute each tile in succession, finally com¬ 
bining the pieces into one large picture. 

There is also scope for video work, moving 
and rotating the viewing plane through the 4D 
parameter space, with different shapes forming 
and collapsing over time. My rough benchmarks 
take 5-10 seconds per frame at 1920 x 1080, so 
for now I’ll wait until faster cheaper graphics 
cards become available. 

6 Obtaining the Implementation 

I wrote the implementation on GNU/Linux De- 
bian Wheezy using proprietary drivers running 
on a quad-core AMD64 processor with NVIDIA 
GTX 550Ti graphics card. 

The source code for the im¬ 
plementation is available at: 

https://gitorious.org/maximus/lyapunov-fm 
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